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Quantum  Mechanical  Calculations  of  the  Second  Virial 
Coefficients  for  Hydrogen* 

Marjorie  E.   Boyd  and  Sigurd  Y.   Larsen 

Quantum  mechanical  second  virial  coefficients  of  dilute 
ortho  and  para  hydrogen  have  been  calculated  for  the  case 
of  a  Lennard- Jones  potential.     These  computations  cover 
the  temperature  range  from  0  °K  to  80  °K.     Special  at- 
tention is  paid  to  the  behaviour  of  the  Boltzmann  and  ex- 
change contributions.     The  purely  statistical  contribution 
to  the  ortho-para  difference  is  shown  to  vanish  to  three 
decimal  places  by  about  10.  5  °K.     Comparison  with  ex- 
perimental data  indicates  that  considerable  improvement 
in  the  potential  function  and  inclusion  of  non-spherical 
terms  will  be  required  at  the  low  temperatures  considered 
here. 

Key  Words  -  Second  Virial  Coefficient,Quantum  Mechanics, 
Hydrogen,   Equation  of  State,   Statistical  Mechanics, 
Intermolecular  Potential 


1.     INTRODUCTION 

This  note  summarizes  the  initial  calculations  in  a  program  aimed  at 
accurate  computation  of  the  quantum  mechanical  second  virial  coef- 
ficient of  dilute  hydrogen.     Ultimately,   it  is  planned  to  carry  out  the 
computations  over  the  complete  range  of  temperatures  from  absolute  zero 
up  to  the  classical  region.     The  work  reported  here  extends  up  to  80  °K. 
Both  the  ortho  and  para  forms  are  considered  as  well  as  their  mixtures. 
While  present  calculations  are  concerned  with  spherical  potential  fields 
for  hydrogen,   it  is  hoped  to  extend  them  eventually  to  include  the  ef- 
fects of  non-spherical  potential  terms. 

One  of  the  aims  of  this  work  is  to  provide  a  sound  theoretical  basis 
for  the  evaluation  and  correlation  of  existing  data  on  hydrogen  now  being 
carried  out  by  the  Cryogenic  Data  Center  at  the  National  Bureau  of 
Standards  Boulder  laboratory.     Our  studies  would  also  aid  in  supplying 
information  in  ranges  where  the  data  are  limited  in  accuracy  or  non- 
existent. 

♦Supported  in  part  by  NBS,   Boulder  Laboratory,    Cryogenic  Data  Center 
under  NASA  Order  R-45,   entitled  "Investigations  and  Studies  of  Liquid 
Hydrogen  and  Other  Cryogenic  Fluids'  ,   for  the  Office  of  Research  Grants 
and  Contracts  of  the  National  Aeronautics  and  Space  Administration. 
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These  calculations  can  be  considered  as  part  of  a  larger  program 
aimed  at  searching  for  an  optimum  intermolecular  potential  function 

suitable  for  determining  not  only  the  second  virial  coefficient  but  also 

the  transport  properties  of  hydrogen.    In  such  an  undertaking  accurate 

calculations  for  the  Lennard- Jones  potential  seem  the  logical  starting 

point  due  to  its  widespread  use. 

2.     BACKGROUND 

The  best  previous  quantum  mechanical  calculation  of  the  second 
virial  coefficient  for  H    is  to  be  found  in  the  work  of  Cohen,  Offerhaus, 
van  Leeuwen,  Roos,  and  de  Boer  '™  (which  also  contains  calculations  of 
the  various  transport  properties),  hereafter  referred  to  as  COLRB.    Their 
work  is  based  on  a  spherically  symmetric  potential  field  assumed  the 
same  for  all  states  of  the  molecule  so  that  the  para-para,  ortho-ortho 
and  para-ortho  interactions  are  all  represented  by  the  same  molecular 
field.     The  second  virial  coefficient  is  then  expressed  in  terms  of  in- 
tegrals which  involve  phase-shifts  arising  from  the  characterization  of 
the  asymptotic  behaviour  of  the  wave  function  of  relative  motion  of  two 
H    molecules.     They  chose  as  the  potential  field  the  well-known 
Lennard-Jones  form 


r  12  6  ) 

4?<r)  =  4ej^         -     (-2-)    J 


where  e  is  the  depth  of  the  potential  well,   and  <j  is  the  "size"  of  the 
molecule  (the  radius  at  which  the  potential  is  zero). 

The  values  used  for  these  parameters  were  obtained  previously  from 
high  temperature  data  by  de  Boer'2)  who  found 

e/K    =    37.00  °K 

N  a   =    1  5. 12  cm    /mol 

from  which  he  obtained,  using  the  1952  values  of  Dumond  and  Cohen  for 
Avogadro's  number  N,  and  the  Boltzmann  constant  K,  the  values 

€     =     31.  88  X  10~4ev 

o 
CT    =       2.928A 

The  Schroedinger  equation  for  the  relative  motion  of  two  particles  in 
this  potential  field  must  then  be  integrated  numerically  for  each  value  of 
energy  (wave  number  k)  and  angular  momentum  (C)  out  to  a  sufficiently 


large  value  of  the  radius  to  determine  the  asymptotic  behaviour  and  thus 
the  phase  shifts  (77    (k)).     The  presence  of  bound  states,   if  any,   and 
their  values  must  be  ascertained.     In  their  calculations  COLRB  used 
values  for  the  phase  shifts  obtained  by  Van  Wijngaarden  and  Duijvestijn 
using  an  electronic  computer  and  calculating  to  three  decimal  places. 
They  obtained  these  phase  shifts  for  values  of  k*,   the  reduced  wave 
number  (l?s  ctv^lE   ,   ^g-        )t   from  0  to  4.  0  at  intervals  of  0.1,   and 
from  4.  0  to  8.  0  at  intervals  of  0.  2  and  for  values  of  the  angular 
momentum  quantum  number, £,   from  0  to  23. 

The  potential  field  they  chose  was  found  to  contain  two  discrete 
energy  levels  and  these  were  determined  to  be  as  follows: 

<t=o     ;    -  3.147  X10"^eV  =  -3.651°K 
1=1:     -  1.227  X  10      eV=  -1.424°K 


For  these  two  values  of  the  angular  momentum  the  phase  shifts  go  to  the 
value  tt  instead  of  zero  as  k*  goes  to  zero. 

The  second  virial  coefficients,   B*  =  B/Ncr    were  computed  for  para- 
hydrogen  (pH2)  and  ortho-hydrogen  (0H2)  in  terms  of  virials  for  Bose  and 
Fermi  statistics,   using  the  following  relations  derived  by  consideration 
of  the  degeneracies  occuring  in  the  two  forms: 


(la)  B*  =       B* 

PH2 


(lb,  B*  =      YBt     +     TBa* 


where  .#3 

ZA 


s        +     '       id 


wi^n  /Vy2  fr  (ZM)r,i  wj  y  iy 


a  tt       T*" 

I  =  even  nos. 
-t  =  odd  nos. 

-372  t*3/2       1_{ZI  +  i)  (e      *  -d 

l  =  even  nos. 
I  =  odd  nos. 
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in  which  the  subscripts  s  or  a  and  the  sums  over  ^-even  or  t-odd  refer 
to  Bose-Einstein  (symmetric)  or  Fermi-Dirac  (antisymmetric)  statistics, 
respectively.     The  various  reduced  quantities  involved  are: 

*  .     . 

A      =     h/CT  Vme 


T       = 

T    (K/c) 

* 

Ei  - 

Vf 

y2   = 

\z  k2  /mKT 

The  first  term,  B       ,   is  the  second  virial  coefficient  of  an  ideal 
quantum  gas  and  has  the  magnitude  A       /  16tt    '      T     3'2.     The  second 
term  is  the  contribution  of  the  continuum  states  of  molecules  interacting 
through  a  potential,   while  the  last  term  gives  the  contribution  of  the 
bound  states  of  energies  E«. 

For  mixtures  of  para  and  ortho  hydrogen  these  authors  use  the 
expression 

B*mix "  *p    V,   +    2*pXo  VoH      +    *o     BoH, 

where  B*        TT    =  1/2  (B       +  B     )  and  the  y's  are  the  percentages  of  each 
p-oH  s  a  A 

species  in  the  mixture. 

Results  are  given  at  temperatures  of  1.  85*K  and  at  intervals  of  3.  7* 
from  3.  7*K  to  37*K.     (See  table  I)    Comparison  with  experiments  in- 
dicated that  their  results  were  about  10%  too  large  at  low  temperatures. 
COLRB  believed  the  discrepancy  to  be  due  largely  to  the  lack  of 
knowledge  of  the  extent  of  influence  of  non-spherical  terms  in  the 
potential  over  the  temperature  range  involved.     To  investigate  this 
possibility  further  these  authors  applied  the  Wang  Chang  '    '  classical 
non-spherical  corrections  to  the  high  temperature  theoretical  curve. 
They  then  redetermined  the  Lennard-Jones  parameters  from  the  re- 
sulting curve.     The  low  temperature  discrepancy  between  theory  and 
experiment  was  then  reduced  to  3%.     However,  they  now  considered 
the  parameters  to  be  of  doubtful  meaning  and  the  procedure  was 
regarded  as  unsatisfactory  since  the  non-spherical  corrections  should 
also  be  applied  at  low  T. 

3.     CALCULATIONS 

Our  initial  calculations  for  the  quantum  mechanical  second 
virial  coefficients  are,   in  general,   an  amplification  of  the  computations 
of  COLRB,   with  several  alterations  and  improvements.    We  employ  a 
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formulation  in  terms  of  the  phase  shifts,  as  they  have  done,   but  have 
used  a  different  arrangement  of  Eq.  (1)  in  which  we  have  separated  the 
contribution  of  a  Boltzmann  gas  (B- direct)  from  the  effects  due  to 
quantum  statistics  (B-exchange).  ^' 

For  a  gas  of  spin  S  we  may  then  write  the  virial  as  a  sum  of  two 
terms: 

direct        exch 


where: 


XT    .  2 

k 


<3'  Bdirect=  -2l/X   PU+1)     7"  /  dk   k    V1"  e 


2    TT 


and 


(4)  B___  =  + 


-  21/2NX3 


*£(V(2«)j^L/,dkk 


exch       (2S+1) 

3 


n  V         2D/^r25 


Here  X^  =    h/  ^2irmKT      ,    (3  =  l/KT,   and  E      ,  is  the  energy  of  the 

th  n' 

n      bound  state   for  angular  momentum  £.     In  B  the  upper  signs 

refer  to  the  symmetrical  (Bose)  statistics  while  the  lower  signs  refer 
to  the  Fermi  or  antisymmetric  case. 

(4) 

Analysis  of  the  behaviour  of  the  B        ,    term         as  well  as 

exch 

numerical  evaluations  for  various  cases' ^>  6)  demonstrates  that  the 

quantum  statistical  effects  embodied  in  it  are  sharply  suppressed  with 

rising  temperature  whenever  a  strong  repulsive  core  is  present  in  the 

potential.     In  fact,   above  a  few  degrees  Kelvin  the  Lennard-Jones  gas 

may  be  regarded  as  effectively  a  Boltzmann  gas. 

This  equation  also  has  the  advantage  of  displaying  the  very 

simple  dependence   of  the  second  virial  coefficient  on  the  spin.     In 

addition,   the  subtraction  involved  between  the  interacting  gas  and  the 

ideal  quantum  gas  contributions  in  evaluating  the  B       ,   term  makes 

exch 


the  effect  of  numerical  inaccuracies  quite  conspicuous.    That  is,  they 
cause  the  subtraction  to  be  inexact  as  B       ,    -»  0  and  erratic  results 

occur  of  the  order  of  magnitude  of  the  digit  in  error. 

Turning  our  attention  now  to  hydrogen  we  first  consider  the  case 
of  the  para  modification.     The  spins  of  the  two  hydrogen  atoms  in  the 
molecule  are  opposed  giving  a  net  spin  of  zero  and  in  its  ground  state 
the  rotational  quantum  number  J  is  zero.    We  thus  have  no  degeneracy 
and  only  symmetric  wave  functions  are  allowed.     That  is,  we  have  a 
spinless  Bose  gas    and 


<5>  B  „      =    B^       f    +    BB°S°    (8=0) 
pH  direct  exch 

where  B  "  ,    is  calculated  setting  S  =  0  and  using  the  upper  signs  in 
_       ...  exch 
Eq.    (4). 

For  ortho  hydrogen  the  situation  is  more  complicated.     Since 
the  spins  of  the  atoms  are  parallel  they  must  be  added,    giving  a  spin 
of  1  for  the  H    molecule.     Thus,  we  have  a  3-fold  degenerate  state, 

where  the  z  component  of  spin,    of  interest  in  a  collision  of  two 

molecules,    can  be  1,    0,    or  -1.     in  addition,  the  ground  state  of  the 

ortho  molecule  has  rotational  quantum  number  J  =  1,   also  three  fold 

degenerate.    The  result  is  thus  9  possible  states  for  the  molecule  and 

9     =  81  possibilities  when  two  molecules  interact.     Considering 

M  =  J  +  S  for  each  molecule  we  have  the  possible  M  values  of  2  (which 

is  5-fold  degenerate),  1  (3-fold  degenerate),   or  0.     Collisions  between 

molecules  of  different  M  values  will  give  only  Boltzmann  type  terms, 

so  no  contributions  will  occur  to  b       ,    in  these  cases.     For  two 

exch  o 

molecules  of  M  =  2  (degeneracy  5)  we  get  bL  -  25  cases,   between 

those  of  M  =  1,   32  =  9  cases,  and  between  those  of  M  =  0,   1  case. 

Thus  for  the  exchange  term,  we  get,  averaging  over  all  the  possible 

cases: 


B  .    RB0Se  ii_        r25e  1 +   qc—i +  _i 1 

exch,  oH2    exch  l°  u;  ^81    L    (2S+l)g=2   ^  (2S+l)g=1  (2S+1), J 


S=0i 


bEosI    (s=o)  :?_ 


exch 
81 


1 


1 


BB°S!         (S=0) 


9  exch 


and  the  total  second  virial  coefficient  for  ortho  hydrogen  will  be 


Bose     ,„   „» 
(6)    B   TT       =      B,  +        1  B        7       S=0 

v   '      oH  direct  — —         exch 

Numerical  improvements  over  the  previous  work  of  COLRB  are  as 
follows.     First,   the  phase  shifts  are  calculated  to  five  figure  accuracy 
in  most  cases  instead  of  the  three  obtained  by  the  previous  workers. 
Second,   instead  of  cutting  off  the  sum  over  I  at  an  arbitrary  value  (23 
in  the  previous  work),   we  have  included  the  phase  shifts  for  all  values 
of  £  .     This  is  done  by  finding  for  each  k*  the  value  of  I  at  which  the 
phase  shift  computed  by  integration  of  the  Schroedinger  equation  agrees 
with  the  value  obtained  by  the  Born  approximation  to  within  one  digit  in 
the  4th  decimal  place.     The  sum  over  the  Born  approximation  values  of 
the  phase  shifts  from  this  value  of  £  to  t=  »can  then  be  evaluated 
analytically. 

The  third  numerical  .improvement  was  to  increase  the  range  of  k* 
values  and  to  decrease  the  interval  size  between  k*  values  used  in  the 
numerical  integration  of  Eqs.   (3)  and  (4).     We  have  calculated  phase 
shifts  up  to  k*  =  16,   twice  the  range  of  COLRB  at  intervals  of  1/16,  in- 
terpolating the  sums  over  them  to  1/64.     This  has  the  effect  of  improving 
the  accuracy  of  the  virial  coefficients  obtained  at  the  low  temperatures 
and  of  allowing  their  calculation  at  higher  temperatures  than  previously 
obtained.    We  give  values  to  80 °K. 

The  values  of  the  energies  of  the  bound  states  in  this  molecule 
were  obtained  using  a  general  computer  program  for  calculating  bound 
states  for  a  Lennard-Jones  potential  of  any  given  parameters  which  was 
written  by  J.  E.  Kilpatrick,   M.   Fraser  and  S.  Y.   Larsen( unpublished). 
Using  the  de  Boer  parameters  these  energies  were  calculated  to  seven 
digits,   as  contrasted  to  the  four  digits  used  by  COLRB.     This  is 
necessary  in  order  to  obtain  more  digits  in  the  low  temperature  values 
of  the  virial  as  can  be  easily  seen  by  examining  the  individual  terms 
given  in  Table  II.     The  lower  bound  state  was  determined  to  be  at 
-3.651296°K,  agreeing  with  the  COLRB  value  to  the  four  digits  they 
used.     The  upper  one  we  have  computed  is  at  -1.422061°K,  which  differs 
from  the  COLRB  value  in  the  fourth  digit.     The  behaviour  of  the  exchange 
term  is  especially  sensitive  to  these  values  and  displays  conspicuously 


the  effect  of  an  error.     In  this  case,  the  erroneous  digit  caused  B 

to  change  sign  when  it  became  of  the  order  of  magnitude  of  the  decimal 

digit  affected  by  the  error. 

Examination  of  the  results  obtained  indicates  general  agreement 

with  the  results  of  COLRB,   as  one  would  expect.     Table  I  gives  our 

results  for  para  and  ortho  H     compared  to  COLRB' s  calculations.    At 

the  lowest  temperature  shown,   1.  85°K,    our  result  for  para  H    is  about 

10  cm  /mole  lower  (only  about  0.1%)  than  the  earlier  calculations,   while 

at  the  upper  end  of  their  range,    37°K,     the  change  is  about  1.  6cm3/mole 

or  about  3%  change  in  the  direction  of  the  experimental  values,   which 

are  some  7%  to  8%  lower  still.     It  is  interesting  to  note  that   7%  is 

about  the  amount  of  improvement  that  was  made  in  the  calculations  of 

COLRB     by-employing  the  classical  non-spherical  correction  at  high  T. 

Thus,    if  applied  to  our  curve,   agreement  with  experiment  within  about 

a  per  cent  or  so  would  result,   but  this  must  be  regarded  as  rather 

fortuitous  at  the  present  stage  of  computation.     Table  III  gives  our 

results  for  B  ,.       ^,    B^ose     and  the  ortho-para  difference  up  to  80 °K. 
direct       exch 

The  part  of  the  ortho-para  difference  which  is  due  to  statistics  alone 

(thst  is,   the  potential  function  is  assumed  the  same  for  both)  depends 

only  on  B        ,    which  can  be  seen  to  have  vanished  to  the  accuracy  we 

exch 
have  calculated  by  about  10.  5°K.     We  emphasize,   however,   that  if  the 

potential  function  does  differ  for  the  different  rotational  states,   as  is 

the  case  in  reality,   the  ortho-para  difference  would  persist  to  higher 

temperatures. 

It  should  be  noted  here  that  at  least  two  re-determinations  of  the 
Lennard-Jones  parameters  for  H    have  been  made  since  the  calculations 
of  COLRB.    Michels,    De  Graaff  and  Ten  SeldamC)   made  a  new  fitting  of 
these  parameters  based  on  experimental  data  obtained  in  their  labora- 
tory'8' for  the  temperature  range  -175°C  to  +150°C.     It  was  not  pos- 
sible to  find  a  set  of  values  for  the  parameters  which  would  enable  the 
fitting  of  the  virials  within  the  experimental  accuracy  over  the  whole 
temperature  range  involved,   even  with  the  introduction  of  high  tempera- 
ture corrections  for  the  non-spherical  effects.     However,   the  best  set 
of  values  was  given  as 

e/K    =     36.7°K 

3  3 

Nor     =15.  60  cm    /  mol 

It  should  be  noted  that  the  expression  for  the  second  virial  coefficient 
at  high  T  used  by  these  authors  is  incorrect,    containing  a  term  for  the 
ideal  gas  (statistical)  contribution  which  can  be  shown  to  have  been 
cancelled  out  by  the  two-body  exchange  effects  at  this  temperature^4). 
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It  is  not  likely  that  this  would  account  for  the  whole  difficulty  in 
obtaining  a  fit  over  the  range  of  temperatures  involved,    but  the  para- 
meter values  given  by  these  authors  must  be  considered  somewhat  un- 
satisfactory. 

(9) 

More  recently,    Diller  and  Mason         fitted  viscosity  data  and 

obtained  the  values 


c/K    =     37.2°k 
a       =      2.97A 

using  mainly  the  date  of  Coremans,   et.  al. ,  and  Johnston  and 

McCloskey  '^'  but  also  measurements  of  van  Itterbeek,   et.   al.,   '     ' 
and  Becker  and  Misenta  U^),     It  has  been  known  for  some  time,    for 
example     in  the  case  of  He,   that  the  best  parameters  or  even  forms  of 
potentials   for  fitting  the  second  virial  are  not  the  best  for  the 
viscosity,  U^)  so  it  is  not  surprising  that  the  values  of  Diller  and  Mason 
are  different  from  those  of  Michels,   et.   al.     The  main  conclusion  about 
them  which  can  be  drawn  at  present  is  that  the  question  of  optimum 
potential  parameters  for  H     is  still  an  open  one.     In  fact,    more 
recently  H.   Hanley  and  M.  Klein  of  NBS  (unpublished  work)  have  ex- 
amined the    possibility  of  fitting  the  data  over  a  large  temperature 
range  in  the  classical  region  by  means  of  various  three -para  meter 
potentials  and  concluded  that  it  was  not  possible  to  do  so  with  any 
of  them. 

In  order  to  extend  our  calculations  over  the  complete  range  of 
temperatures  we  have  recently  computed  new  phase  shifts  extending 
out  to  a  value  of  kcr  =  42.  5.  These  should  allow  the  computation  of 
the  second  virial  coefficient  up  to  the  classical  region,, 


The  computer  programs  used  in  these  calculations  are  an 
adaption  and  elaboration  of  codes  originally  developed  by  Kilpatrick, 
Keller,   Hammel  and  Metropolis  17  who  were  interested  in  helium. 
These  new  codes,  which  we  have  used  for  the  calculations  reported 
in  this  paper,  were  developed  by  the  authors  in  collaboration  with 
Professor  Kilpatrick  and  the  support  of  the  Los  Alamos  Scientific 
Laboratory.    We  are  especially  grateful  to  Malcolm  Fraser, 
Roger  Lazarus  and  Kathleen  Witte  of  LASL  for  much  of  the  programming 
and  the  extensive  testing  of  numerical  methods  as  well  as  for  their 
enthusiastic  aid  and  cooperation. 
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Table  lb 
Second  Virial  Coefficients  of  Hydrogen  -  Boltzmann  Region 


IJI* 

T°K 

B   (cm3 

/mole ) 

C0LRB  Calculation0 

This   Calculation 

Experimental 

.1*0 

lU.8 

-201.1 

-20U.il 

-237     a 

M 

16.65 

-173.17 

-191     a 

.50 

18.5 

-11+6.8 

-11+9.1+2 

-171*     a 

.55 

20.35 

-130.65 

-152     a 

.60 

22.-2 

-113.2 

-115.M 

-127. 8b 

.65 

2k,  05 

-102.886 

-112.5 

.70 

25.9 

-  90.3 

-  92.31*8 

-100.9 

.75 

27.75 

-  83.383 

-  90.90 

.80 

29.6 

-  73.8 

-  75.666 

-  82.1*1 

.85 

31.1+5 

-  68.956 

-  75.03 

.90 

33.3 

-  61.  h 

-  63.070 

-  68.52 

.95 

35.15 

-  57.866 

-  62.72 

1.00 

37.0 

-  51.6 

-  53.233 

-  57.51+ 

1.081 

1*0.0 

-  U6.716 

-  50.32 

1.216 

1+5.0 

-  37.950 

-  1+0.80 

1.351 

50.0 

-  31.08 

-  33.39 

1.1*86 

55.0 

-  25.56 

-  27.1+8° 

1.622 

60.0 

-   21.02 

-   22.70 

1.757 

65.O 

-  17.2 

-  18.6U 

I.892 

70.0 

-  ll+.O 

-  15.22 

2.027 

75.0 

-  11. 

-  12.1*2 

2.162 

80.0 

-     9. 

-     9.88 

a) 
b) 
c) 


Values  interpolated  from  data  in  Reference  15. 
Data  from  this  temperature  up  from  Reference  16. 
Data  from  Reference  1. 
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Table  III 


Direct  and  Exchange  Second  Virial  Coefficients  for  H, 


T(°K) 

direct  (cm  /mole) 

Bose 
gxch 

(B0H2  "  V,' 

1.85 

-871+9.6 

-12U6.9 

+1108.3 

3.70 

-177^.61 

-  18. Ik 

+    16.13 

5.55 

-  869.086 

.870 

+   0.773 

l.k 

-  552.828 

.065 

+   O.O58 

9.25 

-  397. II1* 

.006 

+   0.005 

11.1 

-  305.707 

.000 

0.000 
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